;solution for euler 181


;;References
;http://www2.toki.or.id/book/AlgDesignManual/BOOK/BOOK4/NODE143.HTM
;http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
;http://mitpress.mit.edu/sicp/full-text/sicp/book/node20.html
;http://sandbox.mc.edu/~bennet/cs220/codeex/pgenintro.html
;http://www.jstor.org/pss/2305069

;from euler7
(defun isprime (n lowerprimes)
	(loop for prime in lowerprimes do
		(when (> (* prime prime) n) (return-from isprime t))
		(if (= (mod n prime) 0)
			(return-from isprime nil)))
	(return-from isprime t))

;from euler7
(defun prime-gen (num)
	(let ((primes (list 2)) (i 1) (count 1))
		(loop
			(when (> i num) (return))
			(setf i (+ i 2))
			(when (isprime i primes)
				(incf count)
				(setf (cdr (last primes)) (cons i nil))))
		(return-from prime-gen primes)))

(defparameter *trivial-prime-limit* (floor (sqrt (* 6 (expt 10 7)))))
(defparameter *small-primes* (prime-gen *trivial-prime-limit*)) 
		

(defun nth-prime (n)
	(let ((prime *small-primes*))
		(loop for i from 1 to (- n 1) do
			(setf prime (cdr prime)))
		(return-from nth-prime (car prime))))
		
(defun trivial-prime-count (x)
	(let ((count 0) (prime *small-primes*))
		(loop
			(when (> (car prime) x) (return))
			(setf prime (cdr prime))
			(incf count))
		(return-from trivial-prime-count count)))

(defun phi (m b)
	(when (= b 0) (return-from phi (floor m)))
	(return-from phi
		(- (phi m (- b 1)) (phi (/ m (nth-prime b)) (- b 1)))))
		
(defun prime-count (x)
	(when (< x *trivial-prime-limit*)
		(return-from prime-count (trivial-prime-count x)))
	(let ((count 0) (mlim (floor (expt x (/ 1 3))))
			(nlim (floor (expt x (/ 1 2))))
			m n s)
		(setf m (trivial-prime-count mlim))
		(setf n (trivial-prime-count nlim))
		(setf s (- n m))
		(incf count (+ 
				(phi x m) 
				(* m (+ s 1)) 
				(/ (* s (- s 1)) 2)
				-1))
		(loop for k from 1 to s do
			(decf count (trivial-prime-count (/ x (nth-prime (+ m k))))))
		(return-from prime-count count)))
		
			



		
(defun solve-euler181 ()
	(let (factors 
			(limit (* 5 (expt 10 7))))
		(setf factors (prime-count limit))
		(format t "Euler 181 ~a:" (/ (* factors (+ factors 1)) 2))))



(defun square (a) (* a a))

(defun cong-expt (base ex m)
	(when (= ex 0) (return-from cong-expt 1))
	(if (= 0 (mod ex 2))
		(return-from cong-expt
			(mod (square (cong-expt base (/ ex 2) m)) m))
		(return-from cong-expt
			(mod (* (cong-expt base (- ex 1) m) base) m))))
			
(defun factor-twos (n)
	(let ((count 0))	
		(loop
			(when (= 1 (logand 1 n)) (return))
			(setf n (ash n -1))
			(incf count))
		(return-from factor-twos (cons n count))))




;(defun solve-euler181 ()
;	(let ((count 10001) (n 104743))
;		(loop
;			(incf n 2)
;			(when (= (mod n 10000) 1) (print n))
;			(when (> n 50000) (return))
;			(if (isprime n) (incf count)))
;		(format t "Euler 181: ~a" (/ (* (+ count 1) count) 2))))
	
;;;Notes

;(defparameter *a-values* (list 2 7 61))

;(defun fermat-test (a n) ;test using fermat's little theorem
;	(let ((x 1))
;		(dotimes (i (- n 1)) (setf x (mod (* x a) n)))
;		(return-from fermat-test (= x 1))))

;flawed
;(defun isprime (n)
;	(dolist (p *small-primes*)
;		(when (= 0 (mod n p)) (return-from isprime nil)))
;	(let (d s)
;		(setf d (factor-twos (- n 1)))
;		(setf s (cdr d))
;		(setf d (car d))
;		(dolist (a *a-values*)
;			(when (= 1 (cong-expt a d n)) (return-from isprime t))
;			(dotimes (r (- s 1))
;				(when (= (- n 1) (cong-expt a (* d (ash 1 r)) n))
;					(return-from isprime t)))))
;	(return-from isprime nil))

		
;old approach - failed due to memory requirements		
;(defun solve-euler181 ()
;	(let ( (unmarked 0))
;		(setq sieve (make-array 100))
;		(loop
;			(loop for i from unmarked to (/ (length sieve) (+ unmarked 2)) do ;
;				(setf (elt sieve (* i (+ unmarked 2))) 
;					(+ (elt sieve (* i (+ unmarked 2))) 1)))
;			(loop for i in sieve do
;				(setq unmarked i)
;				(when (= i 0) (return)))
;			(when (= unmarked (- (length sieve) 1)) (return)))
;		(format t "Euler 181: ~a" (count 2 sieve))))

;(defun isprime (n lowerprimes)
;	(loop for prime in lowerprimes do
;		(if (= (mod n prime) 0)
;			(return-from isprime nil)))
;	(return-from isprime t))


;(defun prime-gen (maxn)
;	(let (primes (list 2))
;		(loop for i from 3 to maxn do
;			(if (isprime i primes)
;				(push i primes)))
;		(return-from prime-gen primes)))

		
		
		
		